load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin
  model = (/"ACCESS-CM2","ACCESS-ESM1-5","BCC-CSM2-MR","CNRM-CM6-1","GFDL-ESM4",\
            "HadGEM3-GC31-LL","IPSL-CM6A-LR","MIROC6","MRI-ESM2-0","NorESM2-LM"/);

  expt  = (/"hist-aer","hist-GHG","hist-nat","historical"/)

  do i=0,3
    ts   = new((/10,70/),"float")
    resi = new((/1,70/),"float")
    resi = 0.

    do j=0,9
        fread   = addfile("data/pr_pc1_1951-2020_JJAS_"+expt(i)+"_"+model(j)+".nc","r")
        data    = fread->ts
        ts(j,:) = dim_avg_n_Wrap(data,0)
        data    = data-conform(data,dim_avg_n_Wrap(data,0),1)

        data!1 = "year"
        data!0 = "ens"

        nn = dimsizes(data)

        if (nn(0).gt.1) then
          data  = data*sqrt(nn(0)/(nn(0)-1)); rescale
          resi := array_append_record(resi,data,0)
        end if

        delete(data)
    end do 

    ts   = ts-conform(ts,dim_avg_n_Wrap(ts,1),0)
    resi = resi-conform(resi,dim_avg_n_Wrap(resi,1),0)
    tm   = dim_avg_n_Wrap(ts,0)

    ts_5yr = new(14,"float")
    do m=0,13
        ts_5yr(m) = avg(tm(m*5:m*5+4)); signal 
    end do 

    nr = dimsizes(resi(1:,:))
    resi_5yr = new((/nr(0),14/),"float")
    do m=0,13
        resi_5yr(:,m) = dim_avg_n_Wrap(resi(1:,m*5:m*5+4),1); noise
    end do 
    
    printVarSummary(resi_5yr)

    foutname = "data/ts_signal_noise_5yr_"+expt(i)+".nc"
    system("rm -rf "+foutname)
    fout=addfile(foutname,"c")
    fout->ts=ts_5yr 
    fout->resi=resi_5yr

    foutname = "data/ts_signal_5yr_"+expt(i)+".asc" 
    asciiwrite(foutname,ts_5yr)

    delete([/ts,resi,tm,ts_5yr,resi_5yr,fout/])
    end do 

end 

    

